Intro

Make an evenly spaced UTM prediction grid with all spatially varying covariates for the diet and the biomass data

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(tidylog)
library(sp)
library(raster)
library(devtools)
library(RCurl)
library(sdmTMB)
library(terra)
library(ncdf4)
library(chron)

# Source code for map plots
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5

# Source code for lon lat to utm
source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/lon-lat-utm.R")

theme_set(theme_plot())

Read data and depth-raster

# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
  rename(X = x, Y = y)
#> Rows: 9376 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (3): country, haul_id, ices_rect
#> dbl (20): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)

Make the grid with depth

First make a grid for the biomass data, then subset that based on the extend of the stomach data

x <- d$X
y <- d$Y
z <- chull(x, y)

coords <- cbind(x[z], y[z])

coords <- rbind(coords, coords[1, ])

plot(coords[, 1] ~ coords[, 2]) # plot data


sp_poly <- sp::SpatialPolygons(
  list(sp::Polygons(list(sp::Polygon(coords)), ID = 1))
  )

sp_poly_df <- sp::SpatialPolygonsDataFrame(sp_poly,
                                           data = data.frame(ID = 1)
                                           )
cell_width <- 3

pred_grid <- expand.grid(
  X = seq(min(d$X), max(d$X), cell_width),
  Y = seq(min(d$Y), max(d$Y), cell_width),
  year = unique(d$year)
  )

ggplot(pred_grid %>% filter(year == 2019), aes(X, Y)) +
  geom_point(size = 0.1) +
  theme_void() +
  coord_sf()
#> filter: removed 850,905 rows (96%), 31,515 rows remaining


sp::coordinates(pred_grid) <- c("X", "Y")

inside <- !is.na(sp::over(pred_grid, as(sp_poly_df, "SpatialPolygons")))

pred_grid <- pred_grid[inside, ]

pred_grid <- as.data.frame(pred_grid)

ggplot(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000)) + 
  geom_point(size = 0.001, alpha = 0.5) +
  NULL
#> filter: removed 512,784 rows (96%), 18,992 rows remaining


plot_map +
  geom_point(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000), size = 0.001, alpha = 0.5) +
  NULL
#> filter: removed 512,784 rows (96%), 18,992 rows remaining
#> Warning: Removed 704 rows containing missing values (geom_point).


# Add lat and lon
# Need to go from UTM to lat long for this one...
# https://stackoverflow.com/questions/30018098/how-to-convert-utm-coordinates-to-lat-and-long-in-r
xy <- as.matrix(pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
#> mutate: changed 531,776 values (100%) of 'X' (0 new NA)
#>         changed 531,776 values (100%) of 'Y' (0 new NA)
v <- vect(xy, crs="+proj=utm +zone=33 +datum=WGS84  +units=m")
y <- project(v, "+proj=longlat +datum=WGS84")
lonlat <- geom(y)[, c("x", "y")]

pred_grid$lon <- lonlat[, 1]
pred_grid$lat <- lonlat[, 2]

ggplot(filter(pred_grid, year == 1999), aes(lon, lat)) + geom_point()
#> filter: removed 512,784 rows (96%), 18,992 rows remaining


# Add depth now to remove islands and remaining land
# https://gis.stackexchange.com/questions/411261/read-multiple-layers-raster-from-ncdf-file-using-terra-package
# https://emodnet.ec.europa.eu/geoviewer/
dep_raster <- terra::rast("data/Mean depth natural colour (with land).nc")
class(dep_raster)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"
crs(dep_raster, proj = TRUE)
#> [1] "+proj=longlat +datum=WGS84 +no_defs"

plot(dep_raster)


pred_grid$depth <- terra::extract(dep_raster, pred_grid %>% dplyr::select(lon, lat))$elevation

ggplot(pred_grid, aes(lon, lat, color = depth*-1)) + 
  geom_point()


pred_grid$depth <- pred_grid$depth*-1

pred_grid <- pred_grid %>% drop_na(depth)
#> drop_na: removed 74,340 rows (14%), 457,436 rows remaining

pred_grid %>%
  filter(year == 1999) %>% 
  drop_na(depth) %>% 
  #mutate(water = ifelse(depth < 0.00000001, "N", "Y")) %>% 
  ggplot(aes(X*1000, Y*1000, fill = depth)) + 
  geom_raster() +
  NULL
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> drop_na: no rows removed


plot_map + 
  geom_point(data = pred_grid, aes(X*1000, Y*1000), size = 0.001) + 
  geom_sf()
#> Warning: Removed 19628 rows containing missing values (geom_point).


plot_map + 
  geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = depth), size = 0.001) + 
  geom_sf()
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> Warning: Ignoring unknown parameters: size
#> Warning: Removed 762 rows containing missing values (geom_raster).

Substrate

First add lat and lon based on X and Y (utm)

# # Read raster
# substrate <- raster("data/substrate_tif/BALANCE_SEABED_SEDIMENT.tif")
# 
# # Reproject the raster to fit the UTM pred grid...
# sr <- "+proj=utm +zone=33  +datum=WGS84 +units=m " 
# 
# # Project Raster... This takes some time
# projected_sub_raster <- projectRaster(substrate, crs = sr)
# 
# plot(projected_sub_raster)
# 
# # Now extract the values from the saduria raster to pred_grid
# pred_grid$substrate <- raster::extract(projected_sub_raster, pred_grid %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000))
# 
# unique(pred_grid$substrate)
# 
# # Plot
# ggplot(pred_grid, aes(X, Y, color = substrate)) + 
#   geom_point()
# 
# factor(sort(unique(round(pred_grid$substrate))))
# 
# pred_grid$substrate <- round(pred_grid$substrate)
# 
# pred_grid <- pred_grid %>% mutate(substrate = ifelse(substrate == 1, "bedrock", substrate),
#                                   substrate = ifelse(substrate == 2, "hard-bottom complex", substrate),
#                                   substrate = ifelse(substrate == 3, "sand", substrate),
#                                   substrate = ifelse(substrate == 4, "hard clay", substrate),
#                                   substrate = ifelse(substrate == 5, "mud", substrate))
# 
# pred_grid <- pred_grid %>% drop_na(substrate)
# 
# # I. Bedrock.
# # II. Hard bottom complex, includes patchy hard surfaces and coarse sand (sometimes also clay) to boulders.
# # III. Sand including fine to coarse sand (with gravel exposures).
# # IV. Hard clay sometimes/often/possibly exposed or covered with a thin layer of
# # sand/gravel.
# # V. Mud including gyttja-clay to gyttja-silt.
# 
# # Plot
# ggplot(pred_grid, aes(X, Y, fill = substrate)) + 
#   geom_raster() + 
#   coord_sf()

Oxygen

For these variables, we add quarter to the prediction grid

# Add quarter
pred_grid_1 <- pred_grid %>% mutate(quarter = 1)
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
pred_grid_4 <- pred_grid %>% mutate(quarter = 4)
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA

dat <- bind_rows(pred_grid_1, pred_grid_4)

# Downloaded from here: https://resources.marine.copernicus.eu/?option=com_csw&view=details&product_id=BALTICSEA_REANALYSIS_BIO_003_012
# Extract raster points: https://gisday.wordpress.com/2014/03/24/extract-raster-values-from-points-using-r/comment-page-1/
# https://rpubs.com/boyerag/297592
# https://pjbartlein.github.io/REarthSysSci/netCDF.html#get-a-variable
# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-scobi-monthlymeans_1664182224542.nc")

print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-scobi-monthlymeans_1664182224542.nc (NC_FORMAT_CLASSIC):
#> 
#>      1 variables (excluding dimension variables):
#>         float o2b[longitude,latitude,time]   
#>             long_name: Sea_floor_Dissolved_Oxygen_Concentration
#>             missing_value: NaN
#>             standard_name: mole_concentration_of_dissolved_molecular_oxygen_in_sea_water
#>             units: mmol m-3
#>             _FillValue: NaN
#>             _ChunkSizes: 1
#>              _ChunkSizes: 523
#>              _ChunkSizes: 383
#> 
#>      3 dimensions:
#>         time  Size:336
#>             axis: T
#>             long_name: Validity time
#>             standard_name: time
#>             units: days since 1950-01-01 00:00:00
#>             calendar: gregorian
#>             _ChunkSizes: 512
#>             _CoordinateAxisType: Time
#>             valid_min: 15721.5
#>             valid_max: 25917.5
#>         latitude  Size:523
#>             axis: Y
#>             standard_name: latitude
#>             long_name: latitude
#>             units: degrees_north
#>             _CoordinateAxisType: Lat
#>             valid_min: 48.49169921875
#>             valid_max: 65.8914184570312
#>         longitude  Size:383
#>             standard_name: longitude
#>             long_name: longitude
#>             units: degrees_east
#>             axis: X
#>             _CoordinateAxisType: Lon
#>             valid_min: 9.01375484466553
#>             valid_max: 30.2357654571533
#> 
#>     24 global attributes:
#>         references: http://www.smhi.se
#>         institution: Swedish Meterological and Hydrological Institute
#>         history: See source and creation_date attributees
#>         Conventions: CF-1.5
#>         contact: servicedesk_cmems@mercator-ocean.eu
#>         comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#>         bullentin_type: reanalysis
#>         cmems_product_id: BALTICSEA_REANALYSIS_BIO_003_012
#>         title: CMEMS V4 Reanalysis: SCOBI model 3D fields (monthly means)
#>         FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#>         FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#>         FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#>         FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#>         shallowest_depth: 1.50136542320251
#>         deepest_depth: 711.059204101562
#>         source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#>         file_quality_index: 1
#>         creation_date: 2021-11-09 UTC
#>         bullentin_date: 20201201
#>         start_date: 2020-12-01 UTC
#>         stop_date: 2020-12-01 UTC
#>         start_time: 00:00 UTC
#>         stop_time: 00:00 UTC
#>         _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention

# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 9.013755 9.069310 9.124865 9.180420 9.235975 9.291530

lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 48.49170 48.52503 48.55836 48.59170 48.62503 48.65836

# Get time
time <- ncvar_get(ncin,"time")
time
#>   [1] 15721.5 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0
#>  [10] 15994.5 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0
#>  [19] 16267.5 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5
#>  [28] 16541.0 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5
#>  [37] 16816.5 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0
#>  [46] 17090.5 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0
#>  [55] 17363.5 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5
#>  [64] 17637.0 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5
#>  [73] 17912.5 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0
#>  [82] 18185.5 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0
#>  [91] 18459.5 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5
#> [100] 18733.0 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5
#> [109] 19008.5 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0
#> [118] 19281.5 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0
#> [127] 19554.5 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5
#> [136] 19829.0 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5
#> [145] 20104.5 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0
#> [154] 20377.5 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0
#> [163] 20650.5 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5
#> [172] 20924.0 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5
#> [181] 21199.5 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0
#> [190] 21473.5 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0
#> [199] 21746.5 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5
#> [208] 22020.0 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5
#> [217] 22295.5 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0
#> [226] 22568.5 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0
#> [235] 22842.5 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5
#> [244] 23116.0 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5
#> [253] 23391.5 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0
#> [262] 23664.5 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0
#> [271] 23937.5 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5
#> [280] 24212.0 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5
#> [289] 24487.5 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0
#> [298] 24760.5 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0
#> [307] 25033.5 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5
#> [316] 25307.0 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5
#> [325] 25582.5 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0
#> [334] 25856.5 25887.0 25917.5

tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 336
tunits
#> $hasatt
#> [1] TRUE
#> 
#> $value
#> [1] "days since 1950-01-01 00:00:00"

# Get oxygen
dname <- "o2b"

oxy_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(oxy_array)
#> [1] 383 523 336

# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
oxy_array[oxy_array == fillvalue$value] <- NA

# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#>   [1]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
#>  [26]  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2
#>  [51]  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3
#>  [76]  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4
#> [101]  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5
#> [126]  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6
#> [151]  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7
#> [176]  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8
#> [201]  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9
#> [226] 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10
#> [251] 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11
#> [276] 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12
#> [301]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
#> [326]  2  3  4  5  6  7  8  9 10 11 12

index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)

oxy_q1 <- oxy_array[, , index_keep_q1]
oxy_q4 <- oxy_array[, , index_keep_q4]

months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]

years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]

# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(oxy_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(oxy_q4)[3], by = 3)

# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()

oxy_1 <- c()
oxy_2 <- c()
oxy_3 <- c()
oxy_ave_q1 <- c()

oxy_10 <- c()
oxy_11 <- c()
oxy_12 <- c()
oxy_ave_q4 <- c()

# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want. 

for(i in loop_seq_q1) { # We can use q1 as looping index, doesn't matter!
  
  oxy_1 <- oxy_q1[, , (i)]
  oxy_2 <- oxy_q1[, , (i + 1)]
  oxy_3 <- oxy_q1[, , (i + 2)]
  
  oxy_10 <- oxy_q4[, , (i)]
  oxy_11 <- oxy_q4[, , (i + 1)]
  oxy_12 <- oxy_q4[, , (i + 2)]
  
  oxy_ave_q1 <- (oxy_1 + oxy_2 + oxy_3) / 3
  oxy_ave_q4 <- (oxy_10 + oxy_11 + oxy_12) / 3
    
  list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
  list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
  
  dlist_q1[[list_pos_q1]] <- oxy_ave_q1
  dlist_q4[[list_pos_q4]] <- oxy_ave_q4

}

# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)

# Now I need to make a loop where I extract the raster value for each year...

# Filter years in the pred_grid based on quarter
d_sub_oxy_q1 <- dat %>% filter(quarter == 1)
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
d_sub_oxy_q4 <- dat %>% filter(quarter == 4)
#> filter: removed 457,436 rows (50%), 457,436 rows remaining

# Create data holding object
oxy_data_list_q1 <- list()
oxy_data_list_q4 <- list()

# Create factor year for indexing the list in the loop
d_sub_oxy_q1$year_f <- as.factor(d_sub_oxy_q1$year)
d_sub_oxy_q4$year_f <- as.factor(d_sub_oxy_q4$year)

# Loop through each year and extract raster values for the data points
for(i in sort(unique(d_sub_oxy_q1$year_f))) { # We can use q1 as looping index, doesn't matter!
  
  # Set plot limits
  ymin = 54; ymax = 59; xmin = 12; xmax = 22

  # Subset a year
  oxy_slice_q1 <- dlist_q1[[i]]
  oxy_slice_q4 <- dlist_q4[[i]]
  
  # Create raster for that year (i)
  r_q1 <- raster(t(oxy_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  r_q4 <- raster(t(oxy_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r_q1 <- flip(r_q1, direction = 'y')
  r_q4 <- flip(r_q4, direction = 'y')
  
  plot(r_q1, main = paste(i, "Q1"))
  plot(r_q4, main = paste(i, "Q4"))
  
  # Change projection to UTM (same as pred grid)
  sr <- "+proj=utm +zone=33  +datum=WGS84 +units=m " 
  proj_raster_q1 <- projectRaster(r_q1, crs = sr)
  proj_raster_q4 <- projectRaster(r_q4, crs = sr)
  
  # Filter the same years and select only coordinates
  d_slice_q1 <- d_sub_oxy_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
  d_slice_q4 <- d_sub_oxy_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
  
  # Make into a SpatialPoints object
  data_sp_q1 <- SpatialPoints(d_slice_q1, proj4string = CRS(sr))
  data_sp_q4 <- SpatialPoints(d_slice_q4, proj4string = CRS(sr))
  
  # Extract raster value (oxygen)
  rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
  rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")
  
  # Add in the raster value in the df holding the coordinates for the data
  d_slice_q1$oxy <- rasValue_q1
  d_slice_q4$oxy <- rasValue_q4
  
  # Add in which year
  d_slice_q1$year <- i
  d_slice_q4$year <- i

  # Now the unit of oxygen is mmol/m3. I want it to be ml/L. The original model is in unit ml/L
  # and it's been converted by the data host. Since it was converted without accounting for
  # pressure or temperature, I can simply use the following conversion factor:
  # 1 ml/l = 103/22.391 = 44.661 μmol/l -> 1 ml/l = 0.044661 mmol/l = 44.661 mmol/m^3 -> 0.0223909 ml/l = 1mmol/m^3
  # https://ocean.ices.dk/tools/unitconversion.aspx

  d_slice_q1$oxy <- d_slice_q1$oxy * 0.0223909
  d_slice_q4$oxy <- d_slice_q4$oxy * 0.0223909
    
  # Create a index for the data last where we store all years (because our loop index
  # i is not continuous, we can't use it directly)
  index_q1 <- as.numeric(as.character(d_slice_q1$year))[1] - 1992
  index_q4 <- as.numeric(as.character(d_slice_q4$year))[1] - 1992
  
  # Add each years' data in the list
  oxy_data_list_q1[[index_q1]] <- d_slice_q1
  oxy_data_list_q4[[index_q4]] <- d_slice_q4

}

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)


# Now create a data frame from the list of all annual values
big_dat_oxy_q1 <- dplyr::bind_rows(oxy_data_list_q1)
big_dat_oxy_q4 <- dplyr::bind_rows(oxy_data_list_q4)
big_dat_oxy <- bind_rows(mutate(big_dat_oxy_q1, quarter = 1),
                         mutate(big_dat_oxy_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA

Temperature

# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1664183191233.nc")

print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1664183191233.nc (NC_FORMAT_CLASSIC):
#> 
#>      1 variables (excluding dimension variables):
#>         float bottomT[longitude,latitude,time]   
#>             standard_name: sea_water_potential_temperature_at_sea_floor
#>             units: degrees_C
#>             long_name: Sea floor potential temperature
#>             missing_value: NaN
#>             _FillValue: NaN
#>             _ChunkSizes: 1
#>              _ChunkSizes: 523
#>              _ChunkSizes: 383
#> 
#>      3 dimensions:
#>         time  Size:336
#>             axis: T
#>             long_name: Validity time
#>             standard_name: time
#>             units: days since 1950-01-01 00:00:00
#>             calendar: gregorian
#>             _ChunkSizes: 512
#>             _CoordinateAxisType: Time
#>             valid_min: 15721.5
#>             valid_max: 25917.5
#>         latitude  Size:523
#>             axis: Y
#>             standard_name: latitude
#>             long_name: latitude
#>             units: degrees_north
#>             _CoordinateAxisType: Lat
#>             valid_min: 48.49169921875
#>             valid_max: 65.8914184570312
#>         longitude  Size:383
#>             standard_name: longitude
#>             long_name: longitude
#>             units: degrees_east
#>             axis: X
#>             _CoordinateAxisType: Lon
#>             valid_min: 9.01375484466553
#>             valid_max: 30.2357654571533
#> 
#>     24 global attributes:
#>         references: http://www.smhi.se
#>         institution: Swedish Meterological and Hydrological Institute
#>         history: See source and creation_date attributees
#>         Conventions: CF-1.5
#>         contact: servicedesk_cmems@mercator-ocean.eu
#>         comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#>         bullentin_type: reanalysis
#>         cmems_product_id: BALTICSEA_REANALYSIS_PHY_003_011
#>         title: CMEMS V4 Reanalysis: NEMO model 3D fields (monthly means)
#>         FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#>         FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#>         FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#>         FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#>         shallowest_depth: 1.50136542320251
#>         deepest_depth: 711.059204101562
#>         source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#>         file_quality_index: 1
#>         creation_date: 2021-11-09 UTC
#>         bullentin_date: 20201201
#>         start_date: 2020-12-01 UTC
#>         stop_date: 2020-12-01 UTC
#>         start_time: 00:00 UTC
#>         stop_time: 00:00 UTC
#>         _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention

# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 9.013755 9.069310 9.124865 9.180420 9.235975 9.291530

lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 48.49170 48.52503 48.55836 48.59170 48.62503 48.65836

# Get time
time <- ncvar_get(ncin,"time")
time
#>   [1] 15721.5 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0
#>  [10] 15994.5 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0
#>  [19] 16267.5 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5
#>  [28] 16541.0 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5
#>  [37] 16816.5 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0
#>  [46] 17090.5 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0
#>  [55] 17363.5 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5
#>  [64] 17637.0 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5
#>  [73] 17912.5 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0
#>  [82] 18185.5 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0
#>  [91] 18459.5 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5
#> [100] 18733.0 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5
#> [109] 19008.5 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0
#> [118] 19281.5 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0
#> [127] 19554.5 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5
#> [136] 19829.0 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5
#> [145] 20104.5 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0
#> [154] 20377.5 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0
#> [163] 20650.5 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5
#> [172] 20924.0 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5
#> [181] 21199.5 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0
#> [190] 21473.5 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0
#> [199] 21746.5 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5
#> [208] 22020.0 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5
#> [217] 22295.5 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0
#> [226] 22568.5 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0
#> [235] 22842.5 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5
#> [244] 23116.0 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5
#> [253] 23391.5 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0
#> [262] 23664.5 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0
#> [271] 23937.5 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5
#> [280] 24212.0 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5
#> [289] 24487.5 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0
#> [298] 24760.5 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0
#> [307] 25033.5 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5
#> [316] 25307.0 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5
#> [325] 25582.5 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0
#> [334] 25856.5 25887.0 25917.5

tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 336
tunits
#> $hasatt
#> [1] TRUE
#> 
#> $value
#> [1] "days since 1950-01-01 00:00:00"

# Get temperature
dname <- "bottomT"

temp_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(temp_array)
#> [1] 383 523 336

# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
temp_array[temp_array == fillvalue$value] <- NA

# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#>   [1]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
#>  [26]  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2
#>  [51]  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3
#>  [76]  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4
#> [101]  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5
#> [126]  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6
#> [151]  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7
#> [176]  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8
#> [201]  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9
#> [226] 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10
#> [251] 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11
#> [276] 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12
#> [301]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
#> [326]  2  3  4  5  6  7  8  9 10 11 12

index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)

temp_q1 <- temp_array[, , index_keep_q1]
temp_q4 <- temp_array[, , index_keep_q4]

months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]

years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]

# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(temp_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(temp_q4)[3], by = 3)

# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()

temp_1 <- c()
temp_2 <- c()
temp_3 <- c()
temp_ave_q1 <- c()

temp_10 <- c()
temp_11 <- c()
temp_12 <- c()
temp_ave_q4 <- c()

# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want. 

for(i in loop_seq_q1) {
  
  temp_1 <- temp_q1[, , (i)]
  temp_2 <- temp_q1[, , (i + 1)]
  temp_3 <- temp_q1[, , (i + 2)]
  
  temp_10 <- temp_q4[, , (i)]
  temp_11 <- temp_q4[, , (i + 1)]
  temp_12 <- temp_q4[, , (i + 2)]
  
  temp_ave_q1 <- (temp_1 + temp_2 + temp_3) / 3
  temp_ave_q4 <- (temp_10 + temp_11 + temp_12) / 3
  
  list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
  list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
  
  dlist_q1[[list_pos_q1]] <- temp_ave_q1
  dlist_q4[[list_pos_q4]] <- temp_ave_q4
  
}

# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)

# Now I need to make a loop where I extract the raster value for each year...
# The data is called dat so far in this script

# Filter years in the data frame to only have the years I have temperature for
d_sub_temp_q1 <- dat %>% filter(quarter == 1) %>% filter(year %in% names(dlist_q1)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
d_sub_temp_q4 <- dat %>% filter(quarter == 4) %>% filter(year %in% names(dlist_q4)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed

# Create data holding object
temp_data_list_q1 <- list()
temp_data_list_q4 <- list()

# Create factor year for indexing the list in the loop
d_sub_temp_q1$year_f <- as.factor(d_sub_temp_q1$year)
d_sub_temp_q4$year_f <- as.factor(d_sub_temp_q4$year)

# Loop through each year and extract raster values for the data points
for(i in unique(d_sub_temp_q1$year_f)) {
  
  # Set plot limits
  ymin = 54; ymax = 59; xmin = 12; xmax = 22
  
  # Subset a year
  temp_slice_q1 <- dlist_q1[[i]]
  temp_slice_q4 <- dlist_q4[[i]]
  
  # Create raster for that year (i)
  r_q1 <- raster(t(temp_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  r_q4 <- raster(t(temp_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  
  # Flip...
  r_q1 <- flip(r_q1, direction = 'y')
  r_q4 <- flip(r_q4, direction = 'y')
  
  plot(r_q1, main = paste(i, "Q1"))
  plot(r_q4, main = paste(i, "Q4"))
  
  # Change projection to UTM (same as pred grid)
  sr <- "+proj=utm +zone=33  +datum=WGS84 +units=m " 
  proj_raster_q1 <- projectRaster(r_q1, crs = sr)
  proj_raster_q4 <- projectRaster(r_q4, crs = sr)
  
  # Filter the same years and select only coordinates
  d_slice_q1 <- d_sub_temp_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
  d_slice_q4 <- d_sub_temp_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
  
  # Make into a SpatialPoints object
  data_sp_q1 <- SpatialPoints(d_slice_q1)
  data_sp_q4 <- SpatialPoints(d_slice_q4)
  
  # Extract raster value (temperature)
  rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
  rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")
  
  # Now we want to plot the results of the raster extractions by plotting the
  # data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for pl)
  df_q1 <- as.data.frame(data_sp_q1)
  df_q4 <- as.data.frame(data_sp_q4)
  
  # Add in the raster value in the df holding the coordinates for the data
  d_slice_q1$temp <- rasValue_q1
  d_slice_q4$temp <- rasValue_q4
  
  # Add in which year
  d_slice_q1$year <- i
  d_slice_q4$year <- i
  
  # Create a index for the data last where we store all years (because our loop index
  # i is not continuous, we can't use it directly)
  index_q1 <- as.numeric(d_slice_q1$year)[1] - 1992
  index_q4 <- as.numeric(d_slice_q4$year)[1] - 1992
  
  # Add each years' data in the list
  temp_data_list_q1[[index_q1]] <- d_slice_q1
  temp_data_list_q4[[index_q4]] <- d_slice_q4
  
  }

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)


# Now create a data frame from the list of all annual values
big_dat_temp_q1 <- dplyr::bind_rows(temp_data_list_q1)
big_dat_temp_q4 <- dplyr::bind_rows(temp_data_list_q4)
big_dat_temp <- bind_rows(mutate(big_dat_temp_q1, quarter = 1),
                          mutate(big_dat_temp_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA

Salinity

# https://data.marine.copernicus.eu/product/BALTICSEA_REANALYSIS_PHY_003_011/download?dataset=dataset-reanalysis-nemo-monthlymeans

# Open the netCDF file
ncin <- nc_open("data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1668587452211.nc")

print(ncin)
#> File data/NEMO_Nordic_SCOBI/dataset-reanalysis-nemo-monthlymeans_1668587452211.nc (NC_FORMAT_CLASSIC):
#> 
#>      1 variables (excluding dimension variables):
#>         float sob[longitude,latitude,time]   
#>             long_name: Sea water salinity at sea floor
#>             missing_value: NaN
#>             standard_name: sea_water_salinity
#>             units: 0.001
#>             _FillValue: NaN
#>             _ChunkSizes: 1
#>              _ChunkSizes: 523
#>              _ChunkSizes: 383
#> 
#>      3 dimensions:
#>         time  Size:335
#>             axis: T
#>             long_name: Validity time
#>             standard_name: time
#>             units: days since 1950-01-01 00:00:00
#>             calendar: gregorian
#>             _ChunkSizes: 512
#>             _CoordinateAxisType: Time
#>             valid_min: 15751
#>             valid_max: 25917.5
#>         latitude  Size:187
#>             axis: Y
#>             standard_name: latitude
#>             long_name: latitude
#>             units: degrees_north
#>             _CoordinateAxisType: Lat
#>             valid_min: 53.1249580383301
#>             valid_max: 59.3248596191406
#>         longitude  Size:199
#>             standard_name: longitude
#>             long_name: longitude
#>             units: degrees_east
#>             axis: X
#>             _CoordinateAxisType: Lon
#>             valid_min: 11.1248445510864
#>             valid_max: 22.12473487854
#> 
#>     24 global attributes:
#>         references: http://www.smhi.se
#>         institution: Swedish Meterological and Hydrological Institute
#>         history: See source and creation_date attributees
#>         Conventions: CF-1.5
#>         contact: servicedesk_cmems@mercator-ocean.eu
#>         comment: Provided by SMHI as a Copernicus Marine Environment Monitoring Service production unit
#>         bullentin_type: reanalysis
#>         cmems_product_id: BALTICSEA_REANALYSIS_PHY_003_011
#>         title: CMEMS V4 Reanalysis: NEMO model 3D fields (monthly means)
#>         FROM_ORIGINAL_FILE__easternmost_longitude: 30.2357654571533
#>         FROM_ORIGINAL_FILE__northernmost_latitude: 65.8914184570312
#>         FROM_ORIGINAL_FILE__westernmost_longitude: 9.01375484466553
#>         FROM_ORIGINAL_FILE__southernmost_latitude: 48.49169921875
#>         shallowest_depth: 1.50136542320251
#>         deepest_depth: 711.059204101562
#>         source: SMHI reanalysis run NORDIC-NS2_1d_20201201_20201201
#>         file_quality_index: 1
#>         creation_date: 2021-11-09 UTC
#>         bullentin_date: 20201201
#>         start_date: 2020-12-01 UTC
#>         stop_date: 2020-12-01 UTC
#>         start_time: 00:00 UTC
#>         stop_time: 00:00 UTC
#>         _CoordSysBuilder: ucar.nc2.dataset.conv.CF1Convention

# Get longitude and latitude
lon <- ncvar_get(ncin,"longitude")
nlon <- dim(lon)
head(lon)
#> [1] 11.12484 11.18040 11.23596 11.29151 11.34706 11.40262

lat <- ncvar_get(ncin,"latitude")
nlat <- dim(lat)
head(lat)
#> [1] 53.12496 53.15829 53.19162 53.22496 53.25829 53.29162

# Get time
time <- ncvar_get(ncin,"time")
time
#>   [1] 15751.0 15780.5 15811.0 15841.5 15872.0 15902.5 15933.5 15964.0 15994.5
#>  [10] 16025.0 16055.5 16086.5 16116.0 16145.5 16176.0 16206.5 16237.0 16267.5
#>  [19] 16298.5 16329.0 16359.5 16390.0 16420.5 16451.5 16481.0 16510.5 16541.0
#>  [28] 16571.5 16602.0 16632.5 16663.5 16694.0 16724.5 16755.0 16785.5 16816.5
#>  [37] 16846.5 16876.5 16907.0 16937.5 16968.0 16998.5 17029.5 17060.0 17090.5
#>  [46] 17121.0 17151.5 17182.5 17212.0 17241.5 17272.0 17302.5 17333.0 17363.5
#>  [55] 17394.5 17425.0 17455.5 17486.0 17516.5 17547.5 17577.0 17606.5 17637.0
#>  [64] 17667.5 17698.0 17728.5 17759.5 17790.0 17820.5 17851.0 17881.5 17912.5
#>  [73] 17942.0 17971.5 18002.0 18032.5 18063.0 18093.5 18124.5 18155.0 18185.5
#>  [82] 18216.0 18246.5 18277.5 18307.5 18337.5 18368.0 18398.5 18429.0 18459.5
#>  [91] 18490.5 18521.0 18551.5 18582.0 18612.5 18643.5 18673.0 18702.5 18733.0
#> [100] 18763.5 18794.0 18824.5 18855.5 18886.0 18916.5 18947.0 18977.5 19008.5
#> [109] 19038.0 19067.5 19098.0 19128.5 19159.0 19189.5 19220.5 19251.0 19281.5
#> [118] 19312.0 19342.5 19373.5 19403.0 19432.5 19463.0 19493.5 19524.0 19554.5
#> [127] 19585.5 19616.0 19646.5 19677.0 19707.5 19738.5 19768.5 19798.5 19829.0
#> [136] 19859.5 19890.0 19920.5 19951.5 19982.0 20012.5 20043.0 20073.5 20104.5
#> [145] 20134.0 20163.5 20194.0 20224.5 20255.0 20285.5 20316.5 20347.0 20377.5
#> [154] 20408.0 20438.5 20469.5 20499.0 20528.5 20559.0 20589.5 20620.0 20650.5
#> [163] 20681.5 20712.0 20742.5 20773.0 20803.5 20834.5 20864.0 20893.5 20924.0
#> [172] 20954.5 20985.0 21015.5 21046.5 21077.0 21107.5 21138.0 21168.5 21199.5
#> [181] 21229.5 21259.5 21290.0 21320.5 21351.0 21381.5 21412.5 21443.0 21473.5
#> [190] 21504.0 21534.5 21565.5 21595.0 21624.5 21655.0 21685.5 21716.0 21746.5
#> [199] 21777.5 21808.0 21838.5 21869.0 21899.5 21930.5 21960.0 21989.5 22020.0
#> [208] 22050.5 22081.0 22111.5 22142.5 22173.0 22203.5 22234.0 22264.5 22295.5
#> [217] 22325.0 22354.5 22385.0 22415.5 22446.0 22476.5 22507.5 22538.0 22568.5
#> [226] 22599.0 22629.5 22660.5 22690.5 22720.5 22751.0 22781.5 22812.0 22842.5
#> [235] 22873.5 22904.0 22934.5 22965.0 22995.5 23026.5 23056.0 23085.5 23116.0
#> [244] 23146.5 23177.0 23207.5 23238.5 23269.0 23299.5 23330.0 23360.5 23391.5
#> [253] 23421.0 23450.5 23481.0 23511.5 23542.0 23572.5 23603.5 23634.0 23664.5
#> [262] 23695.0 23725.5 23756.5 23786.0 23815.5 23846.0 23876.5 23907.0 23937.5
#> [271] 23968.5 23999.0 24029.5 24060.0 24090.5 24121.5 24151.5 24181.5 24212.0
#> [280] 24242.5 24273.0 24303.5 24334.5 24365.0 24395.5 24426.0 24456.5 24487.5
#> [289] 24517.0 24546.5 24577.0 24607.5 24638.0 24668.5 24699.5 24730.0 24760.5
#> [298] 24791.0 24821.5 24852.5 24882.0 24911.5 24942.0 24972.5 25003.0 25033.5
#> [307] 25064.5 25095.0 25125.5 25156.0 25186.5 25217.5 25247.0 25276.5 25307.0
#> [316] 25337.5 25368.0 25398.5 25429.5 25460.0 25490.5 25521.0 25551.5 25582.5
#> [325] 25612.5 25642.5 25673.0 25703.5 25734.0 25764.5 25795.5 25826.0 25856.5
#> [334] 25887.0 25917.5

tunits <- ncatt_get(ncin,"time","units")
nt <- dim(time)
nt
#> [1] 335
tunits
#> $hasatt
#> [1] TRUE
#> 
#> $value
#> [1] "days since 1950-01-01 00:00:00"

# Get Salinity
dname <- "sob"

sal_array <- ncvar_get(ncin,dname)
dlname <- ncatt_get(ncin,dname,"long_name")
dunits <- ncatt_get(ncin,dname,"units")
fillvalue <- ncatt_get(ncin,dname,"_FillValue")
dim(sal_array)
#> [1] 199 187 335

# Get global attributes
title <- ncatt_get(ncin,0,"title")
institution <- ncatt_get(ncin,0,"institution")
datasource <- ncatt_get(ncin,0,"source")
references <- ncatt_get(ncin,0,"references")
history <- ncatt_get(ncin,0,"history")
Conventions <- ncatt_get(ncin,0,"Conventions")

# Convert time: split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth <- as.integer(unlist(tdstr)[2])
tday <- as.integer(unlist(tdstr)[3])
tyear <- as.integer(unlist(tdstr)[1])

# Here I deviate from the guide a little bit. Save this info:
dates <- chron(time, origin = c(tmonth, tday, tyear))

# Crop the date variable
months <- as.numeric(substr(dates, 2, 3))
years <- as.numeric(substr(dates, 8, 9))
years <- ifelse(years > 90, 1900 + years, 2000 + years)

# Replace netCDF fill values with NA's
sal_array[sal_array == fillvalue$value] <- NA

# Next, we need to work with the months that correspond to the quarters that we use.
# loop through each time step, and if it is a good month save it as a raster.
# First get the index of months that correspond to Q4
months
#>   [1]  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2
#>  [26]  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3
#>  [51]  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4
#>  [76]  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5
#> [101]  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6
#> [126]  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7
#> [151]  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8
#> [176]  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9
#> [201] 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10
#> [226] 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11
#> [251] 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12
#> [276]  1  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1
#> [301]  2  3  4  5  6  7  8  9 10 11 12  1  2  3  4  5  6  7  8  9 10 11 12  1  2
#> [326]  3  4  5  6  7  8  9 10 11 12

index_keep_q1 <- which(months < 4)
index_keep_q4 <- which(months > 9)

sal_q1 <- sal_array[, , index_keep_q1]
sal_q4 <- sal_array[, , index_keep_q4]

months_keep_q1 <- months[index_keep_q1]
months_keep_q4 <- months[index_keep_q4]

years_keep_q1 <- years[index_keep_q1]
years_keep_q4 <- years[index_keep_q4]

# Now we have an array with data for that quarter
# We need to now calculate the average within a year.
# Get a sequence that takes every third value between 1: number of months (length)
loop_seq_q1 <- seq(1, dim(sal_q1)[3], by = 3)
loop_seq_q4 <- seq(1, dim(sal_q4)[3], by = 3)

# Create objects that will hold data
dlist_q1 <- list()
dlist_q4 <- list()

sal_1 <- c()
sal_2 <- c()
sal_3 <- c()
sal_ave_q1 <- c()

sal_10 <- c()
sal_11 <- c()
sal_12 <- c()
sal_ave_q4 <- c()

# Now average by quarter. The vector loop_seq_q1 is 1, 4, 7 etc. So first i is 1, 2, 3,
# which is the index we want.

dim(sal_q1)
#> [1] 199 187  83
dim(sal_q4)
#> [1] 199 187  84

# Hmm, we didn't get the first month in the salinity series... repeat month 2 and fill in so the dimensions are correct
sal_q1 <- sal_q1[,,c(1, 1:83)]

dim(sal_q1)
#> [1] 199 187  84

for(i in loop_seq_q1) {

  sal_1 <- sal_q1[, , (i)]
  sal_2 <- sal_q1[, , (i + 1)]
  sal_3 <- sal_q1[, , (i + 2)]

  sal_10 <- sal_q4[, , (i)]
  sal_11 <- sal_q4[, , (i + 1)]
  sal_12 <- sal_q4[, , (i + 2)]

  sal_ave_q1 <- (sal_1 + sal_2 + sal_3) / 3
  sal_ave_q4 <- (sal_10 + sal_11 + sal_12) / 3

  list_pos_q1 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)
  list_pos_q4 <- ((i/3) - (1/3)) + 1 # to get index 1:n(years)

  dlist_q1[[list_pos_q1]] <- sal_ave_q1
  dlist_q4[[list_pos_q4]] <- sal_ave_q4

}

# Now name the lists with the year:
names(dlist_q1) <- unique(years_keep_q1)
names(dlist_q4) <- unique(years_keep_q4)

# Now I need to make a loop where I extract the raster value for each year...
# The cpue data is called dat so far in this script

# Filter years in the cpue data frame to only have the years I have salinity for
d_sub_sal_q1 <- dat %>% filter(quarter == 1) %>% filter(year %in% names(dlist_q1)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed
d_sub_sal_q4 <- dat %>% filter(quarter == 4) %>% filter(year %in% names(dlist_q4)) %>% droplevels()
#> filter: removed 457,436 rows (50%), 457,436 rows remaining
#> filter: no rows removed

# Create data holding object
sal_data_list_q1 <- list()
sal_data_list_q4 <- list()

# Create factor year for indexing the list in the loop
d_sub_sal_q1$year_f <- as.factor(d_sub_sal_q1$year)
d_sub_sal_q4$year_f <- as.factor(d_sub_sal_q4$year)

# Loop through each year and extract raster values for the cpue data points
for(i in unique(d_sub_sal_q1$year_f)) {

  # Set plot limits
  ymin = 54; ymax = 58; xmin = 12; xmax = 22

  # Subset a year
  sal_slice_q1 <- dlist_q1[[i]]
  sal_slice_q4 <- dlist_q4[[i]]

  # Create raster for that year (i)
  r_q1 <- raster(t(sal_slice_q1), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))
  r_q4 <- raster(t(sal_slice_q4), xmn = min(lon), xmx = max(lon), ymn = min(lat), ymx = max(lat),
                 crs = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

  # Flip...
  r_q1 <- flip(r_q1, direction = 'y')
  r_q4 <- flip(r_q4, direction = 'y')

  plot(r_q1, main = paste(i, "Q1"))
  plot(r_q4, main = paste(i, "Q4"))

  # Change projection to UTM (same as pred grid)
  sr <- "+proj=utm +zone=33  +datum=WGS84 +units=m " 
  proj_raster_q1 <- projectRaster(r_q1, crs = sr)
  proj_raster_q4 <- projectRaster(r_q4, crs = sr)

  # Filter the same year (i) in the cpue data and select only coordinates
  d_slice_q1 <- d_sub_sal_q1 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)
  d_slice_q4 <- d_sub_sal_q4 %>% filter(year_f == i) %>% dplyr::select(X, Y) %>% mutate(X = X*1000, Y = Y*1000)

  # Make into a SpatialPoints object
  data_sp_q1 <- SpatialPoints(d_slice_q1)
  data_sp_q4 <- SpatialPoints(d_slice_q4)

  # Extract raster value (salinity)
  rasValue_q1 <- raster::extract(proj_raster_q1, data_sp_q1, method = "bilinear")
  rasValue_q4 <- raster::extract(proj_raster_q4, data_sp_q4, method = "bilinear")

  # Now we want to plot the results of the raster extractions by plotting the cpue
  # data points over a raster and saving it for each year.
  # Make the SpatialPoints object into a raster again (for pl)
  df_q1 <- as.data.frame(data_sp_q1)
  df_q4 <- as.data.frame(data_sp_q4)

  # Add in the raster value in the df holding the coordinates for the cpue data
  d_slice_q1$sal <- rasValue_q1
  d_slice_q4$sal <- rasValue_q4

  # Add in which year
  d_slice_q1$year <- i
  d_slice_q4$year <- i

  # Create a index for the data last where we store all years (because our loop index
  # i is not continuous, we can't use it directly)
  index_q1 <- as.numeric(d_slice_q1$year)[1] - 1992
  index_q4 <- as.numeric(d_slice_q4$year)[1] - 1992

  # Add each years' data in the list
  sal_data_list_q1[[index_q1]] <- d_slice_q1
  sal_data_list_q4[[index_q4]] <- d_slice_q4

}

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)

#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)
#> filter: removed 441,099 rows (96%), 16,337 rows remaining
#> mutate: changed 16,337 values (100%) of 'X' (0 new NA)
#>         changed 16,337 values (100%) of 'Y' (0 new NA)


# Now create a data frame from the list of all annual values
big_dat_sal_q1 <- dplyr::bind_rows(sal_data_list_q1)
big_dat_sal_q4 <- dplyr::bind_rows(sal_data_list_q4)
big_dat_sal <- bind_rows(mutate(big_dat_sal_q1, quarter = 1),
                          mutate(big_dat_sal_q4, quarter = 4))
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
#> mutate: new variable 'quarter' (double) with one unique value and 0% NA
big_dat_oxy <- big_dat_oxy %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
big_dat_temp <- big_dat_temp %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
big_dat_sal <- big_dat_sal %>% mutate(id_env = paste(year, quarter, X, Y, sep = "_"))
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA

env_dat <- left_join(big_dat_oxy, big_dat_temp) %>%
  left_join(big_dat_sal) %>% 
  dplyr::select(id_env, oxy, temp, sal)
#> Joining, by = c("X", "Y", "year", "quarter", "id_env")
#> left_join: added one column (temp)
#>            > rows only in x         0
#>            > rows only in y  (      0)
#>            > matched rows     914,872
#>            >                 =========
#>            > rows total       914,872
#> Joining, by = c("X", "Y", "year", "quarter", "id_env")
#> left_join: added one column (sal)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 914,872
#> > =========
#> > rows total 914,872
  
dat <- dat %>%
  mutate(id_env = paste(year, quarter, X*1000, Y*1000, sep = "_")) %>% 
  left_join(env_dat) %>%
  dplyr::select(-id_env)
#> mutate: new variable 'id_env' (character) with 914,872 unique values and 0% NA
#> Joining, by = "id_env"left_join: added 3 columns (oxy, temp, sal)
#>            > rows only in x         0
#>            > rows only in y  (      0)
#>            > matched rows     914,872
#>            >                 =========
#>            > rows total       914,872

dat <- dat %>% drop_na(temp) %>% drop_na(oxy) %>% drop_na(sal)
#> drop_na: removed 5,712 rows (1%), 909,160 rows remaining
#> drop_na: no rows removed
#> drop_na: removed 1,008 rows (<1%), 908,152 rows remaining

# Temperature
plot_map_fc + 
  geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = temp)) + 
  facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).


# Oxygen
plot_map_fc + 
  geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = oxy)) +
  facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).


# Salinity
plot_map_fc +
  geom_raster(data = filter(dat, quarter == 4), aes(X*1000, Y*1000, fill = sal)) +
  facet_wrap(~year)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
#> Warning: Removed 21336 rows containing missing values (geom_raster).


pred_grid <- dat

Add ICES areas

# https://stackoverflow.com/questions/34272309/extract-shapefile-value-to-point-with-r
# https://gis.ices.dk/sf/
shape <- shapefile("data/ICES_StatRec_mapto_ICES_Areas/StatRec_map_Areas_Full_20170124.shp")
head(shape)
#>   OBJECTID ID ICESNAME SOUTH WEST NORTH EAST AREA_KM2 stat_x stat_y Area_27
#> 0        1 47     47A0  59.0  -44  59.5  -43     3178  -43.5  59.25  14.b.2
#> 1        2 48     48A0  59.5  -44  60.0  -43     3132  -43.5  59.75  14.b.2
#> 2        3 49     49A0  60.0  -44  60.5  -43     3085  -43.5  60.25  14.b.2
#> 3        4 50     50A0  60.5  -44  61.0  -43     3038  -43.5  60.75  14.b.2
#> 4        5 51     51A0  61.0  -44  61.5  -43     2991  -43.5  61.25  14.b.2
#> 5        6 52     52A0  61.5  -44  62.0  -43     2944  -43.5  61.75  14.b.2
#>           Perc      MaxPer RNDMaxPer AreasList Shape_Leng Shape_Area
#> 0 100.00000000 100.0000000       100    14.b.2          3        0.5
#> 1  84.12674145  84.1267414        84    14.b.2          3        0.5
#> 2  24.99803694  24.9980369        25    14.b.2          3        0.5
#> 3  11.97744244  11.9774424        12    14.b.2          3        0.5
#> 4   3.89717625   3.8971762         4    14.b.2          3        0.5
#> 5   0.28578428   0.2857843         0    14.b.2          3        0.5

pts <- SpatialPoints(cbind(pred_grid$lon, pred_grid$lat), 
                     proj4string = CRS(proj4string(shape)))
#> Warning in proj4string(shape): CRS object has comment, which is lost in output

pred_grid$subdiv <- over(pts, shape)$Area_27

# Rename subdivisions to the more common names and do some more filtering (by sub div and area)
sort(unique(pred_grid$subdiv))
#> [1] "3.d.24"   "3.d.25"   "3.d.26"   "3.d.27"   "3.d.28.2"

pred_grid <- pred_grid %>% 
  mutate(sub_div = factor(subdiv),
         sub_div = fct_recode(subdiv,
                              "24" = "3.d.24",
                              "25" = "3.d.25",
                              "26" = "3.d.26",
                              "27" = "3.d.27",
                              "28" = "3.d.28.1",
                              "28" = "3.d.28.2",
                              "29" = "3.d.29"),
         sub_div = as.character(sub_div)) %>% 
  filter(sub_div %in% c("24", "25", "26", "27", "28", 2)) %>% 
  filter(lat > 54 & lat < 59 & lon < 22)
#> Warning: Unknown levels in `f`: 3.d.28.1, 3.d.29
#> mutate: new variable 'sub_div' (character) with 5 unique values and 0% NA
#> filter: no rows removed
#> filter: no rows removed

# Add ICES rectangles
pred_grid$ices_rect <- mapplots::ices.rect2(lon = pred_grid$lon, lat = pred_grid$lat)

plot_map +
  geom_raster(data = filter(pred_grid, year == 1999), aes(X*1000, Y*1000, fill = oxy)) +
  facet_wrap(~sub_div)
#> filter: removed 875,718 rows (96%), 32,434 rows remaining
#> Warning: Removed 1524 rows containing missing values (geom_raster).


pred_grid <- pred_grid %>% dplyr::select(-subdiv)

Save

# Remove variables and save
pred_grid_93_06 <- pred_grid %>% filter(year < 2007)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining
pred_grid_07_19 <- pred_grid %>% filter(year > 2006)
#> filter: removed 454,076 rows (50%), 454,076 rows remaining

write.csv(pred_grid_93_06, file = "data/clean/pred_grid_(1_2).csv", row.names = FALSE)
write.csv(pred_grid_07_19, file = "data/clean/pred_grid_(2_2).csv", row.names = FALSE)